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In this paper, we predict the existence of low-frequency nonlocal plasmon excitations at the 
vacuum-surface interface of a superlattice of N graphene layers interacting with a thick conducting 
substrate. This is different from graphite which allows inter-layer hopping. A dispersion function is 
derived which incorporates the polarization function of the graphene monolayers (MLGs) and the 
dispersion function of a semi-infinite electron liquid at whose surface the electrons scatter specularly. 
We find that this surface plasmon-polariton is not damped by the particle-hole excitations (PHEs) 
or the bulk modes and separates below the continuum mini-band of bulk plasmon modes. For a 
conducting substrate with surface plasmon frequency ujs — (x}pl\/2, the surface plasmon frequency 
of the hybrid structure always lies below uJs . The intensity of this mode depends on the distance of 
the graphene layers from the surface of the conductor, the energy band gap between the valence and 
conduction bands of MLG and, most importantly, on the number of two-dimensional (2D) layers. 
Furthermore, the hybrid structure has no surface plasmon for a sufficiently large number [N ~ 7) of 
layers. The existence of two plasmons with different dispersion relations indicates that quasiparticles 
with different group velocity may coexist for various ranges of wavelength which is determined by 
the number of layers in the superlattice. 

PACS numbers: 73.21.-b, 03.67.Lx, 71.70.Ej,73.20.Mf, 71.45.Gm, 71.10.C 


I. INTRODUCTION 

Several years ago, there was considerable activity in the study of the electronic properties of layered semiconductor 
superlattices CHS] leading to intriguing transport and optical properties which are the result of quantum mechanical 
effects on the nanoscale. Specifically, it was desirable to learn how a periodic or quasi-periodic array [6] of two- 
dimensional electron liquid (2DEL) layers, would lead to a modification of the response of the charge-density excitations 
to an electromagnetic field. Recently, there has been some investigation regarding the Goulomb excitations of an 
Y-layered superlattice of free-standing monolayer graphene layers (MGLs) [7]. However, what was neglected in 
that study is the interaction between epitaxial layers of graphene and a conducting substrate which may give rise 
to composite plasmon-plasmon resonances. A low-frequency mode emerges below the continuum of modes in the 
limit of a large number of layers and is associated with the air-exposed surface of the graphene-conductor substrate 
combination. As a matter of fact, several authors IEHI31 demonstrated how the electronic response properties of 
graphene-metal composites of Ru and Ni, for example, are much different from free-standing structures. Furthermore, 
these complex carbon/metal interfaces are interesting because of the unusual and fundamental physics regarding their 
electronic and magnetic properties at the 2D interfaces. Examples of these systems occur in intercalated graphite 
na, incommensurate transition metal/graphene [15] and carbide/graphene interfaces [16] . Possible motivation for 
pursuing this area of research is the tunability of graphene plasmons by a substrate which is a promising emerging 
field of graphene-based plasmonics. 

Here, we demonstrate how the intensity of the response of the surface plasmon arising when N graphene layers 
interact with a conducting substrate may be adjusted by changing the layer-substrate separation, the energy band 
gap or the number of 2D layers. It has been shown that such an energy gap may be produced by circularly polarized 
light [niiiHi. We employ our model to determine the plasmon excitation dispersion relation for a single layer on a 
conducting substrate and show how our results may simulate the experimental data recently reported by Politano, 
et al. [snsj. Additionally, this paper investigates multi-layer graphene generally, showing how its new surface mode 
depends on the in-plane wave vector and the critical wave vector qc which marks the onset of damping by bulk modes 
in the miniband continuum. 

The model we use for a superlattice consists of N 2D graphene layers whose planes are perpendicular to the z-axis 
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at 2 ; = ai,a 2 , • • • and the substrate occupies the half-space z < 0. Each graphene layer will be described by an 
energy band structure for Dirac fermions and may be intrinsic, gapped or doped. This model does not allow inter-layer 
hopping and hence it differs from the one we employed to describe the anisotropy of 7r-plasmon dispersion in AA- 
stacked graphite [19] . The screening of an externally applied frequency-dependent potential by the polarized medium 
requires a knowledge of the dielectric function of the structure which we obtain in the random-phase-approximation 
(RPA). We present our method of calculation in Sec. [n| Section [111] is devoted to numerical results and discussion 
of our plasmon dispersion when several graphene layers interact with a conducting substrate. The simulated data 
show how the intensity of the modes depends on the number of layers which are stacked, their distance from the 
conducting surface as well as the energy gap. For N = 1 gapless MLG, when the conductor surface-MLG separation 
exceeds a critical distance dc, the intensity of the surface plasmon in the long wavelength regime is sufficiently high 
to be observable up to some cut-off wave vector Qc of the surface plasmon frequency. Beyond Qc, the intensity of ujc 
is very weak until the plasmon wave vector reaches some value q'^. However, when the surface-MLG separation is less 
that dc, the surface plasmon intensity in the long wavelength regime is weak and the mode only appears at shorter 
wavelengths when the in-plane wave vector gy > q'^. Interestingly, for gapped graphene, uJc is completely suppressed 
when the surface- layer separation is less than dc- We then make some concluding remarks in Sec. |IV| regarding the 
inspiration for our work, a summary and significance of our findings, and what new theoretical formalism is presented 
in our paper. 


II. GENERAL FORMULATION OF THE PROBLEM 

In our formalism, we consider a nano-scale system consisting of an arbitrary number of 2D layers and a thick 
conducting material. The layer may be monolayer graphene or a 2DEL such as a semiconductor inversion layer or 
HEMT (high electron mobility transistor). The graphene layer may have a gap, thereby extending the flexibility of 
the composite system that also incorporates a thick layer of dielectric material. The excitation spectra of allowable 
modes will be determined from a knowledge of the non-local dielectric function e(r,r';co’) which depends on position 
coordinates r, r' and frequency uj or its inverse satisfying J dr' A(r, r'; Cc;)e(r', r"; cj) = d(r, r"). 

In operator notation, the dielectric function for the N 2D layers and a semi-infinite structure is given by [20| 

N N N 

e = i -h ^si + ^2D — ^SI + ^2D ~ ^SI ^ ^2D ’ (^) 

i=l i=l i=l 

where e = K~^ with K the inverse dielectric function satisfying 

N 

K = ksi-Ksi-J24%-K . ( 2 ) 

i=l 

In integral form, after Fourier transforming parallel to the x^-plane and suppressing the in-plane wave number g| | and 
frequency cj, we obtain 


roo POO 

K{zi,Z2) = Ksiizi,Z2) K{zuz')kil{z',z")K{z",Z2) . (3) 

J —00 J —00 

Here, the polarization function for the 2D structure is given by 

/ c» 

dz'" v{z',z"')D^^\z"',z") , (4) 

-00 

where the 2D response function obeys 

D^^\z'",z") = nfl^,{q^\,oj)5iz'" - ai)5{z" - a) (5) 

with 112*^ w) the single-particle in-plane response [2TJ|23l|231126] . The 2D RPA ring diagram polarization function 
for graphene with a gap A may be expressed as 
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1 + ss‘ 


s,s'=± \ 

f{s€k) - /(S'ek+q||) 
f ek|| - s'ek+q,, -fuj-iO+ 


,k- (k + q||) + A^ ' 

^|k+q||| 


( 6 ) 


where g = ^ takes account of valley and spin degeneracy. At T = 0, the Fermi-Dirac distribution function is reduced 
to the Heaviside step function /(e,/i;T ^ 0) = 77+(/i — e), In the long wavelength limit, the real part of is 

given by 


Re nSfeii.o,) = - + 

f 3 /i jih? 1 / 2 /i — 

Atv \ 4/i^ — ^ 2huj^ \2jii hw) j 


and the imaginary part by 


Im n(°^(g||,w) = ^ ri+{2ii - tiLo) . (8) 

Upon substituting this form of the polarization function for the monolayer into the integral equation for the inverse 
dielectric function, we have 


/ oo 

dz' Ksi{zi,z')v{z' - ai)K{ai,Z 2 ) ■ 

-oo 

We now set zi = in Eq. and obtain 

^1%/+n^'^;^.(g||,w) J dz' Ksi{aj,,z')v{z'- aj)^ K{aj, Z 2 ) = Ksi{aj,, Z 2 ) . 

These linear algebraic equations may be solved simultaneously and their solutions expressed in matrix form as 
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( 10 ) 


/K{ai,Z2) 

K{a2,Z2) 

\K{a2,N2) 





Ksi{ai,Z2) \ 
Ksi{a2, Z2) 

KsiidN, Z 2 ) 

/ 


where = DetA1^'^^(g|| , uj) with matrix elements given by 

/ oo 

dz' Ksi{aj,,z')v{z' -aj) 

-OO 

and is the transpose cofactor matrix of w). 

In mean-field theory, we have |25j 


( 11 ) 


( 12 ) 
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FIG. 1: (Color online) Plasmon dispersion relation for a semi-infinite conductor which is Coulomb-coupled to monolayer 
graphene for various surface-to-layer separations. In the panels (a), (c) and (e) on the left-hand side, we present density plots 
of the inverse dispersion function 1 /+ zO^) with peaks corresponding to the undamped plasmon modes. The right 
panels (5), (d) and (/) show numerical solutions for the plasmon branches, both Landau damped and undamped. The distances 
chosen are akr = 0.2, 1.0 and 0.5, correspondingly. All plots are provided for extrinsic graphene (doped) with zero energy 
bandgap. 


where 


Ksi 

= V+(z)ls(z-z') - 


Oil) 


1 + e (g||) 
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+ (q||,3 + 3';w) + z - z';w)] )?+(-3')| - 


(13) 


Oil) ^ = 


(14) 




























5 



0.0 0.6 1.2 a /I 2.4 3.0 

%/kF 





FIG. 2: (Color online) Plasmon excitation spectra for a semi-infinite conductor interacting through the Coulomb interaction 
with N — 2 monolayer graphene sheets located at chosen distances from the surface. The left-hand panels (a), (c) and (e) give 
density plots of the inverse dispersion function 1 /{q\\, u iO'^) with peaks corresponding to the plasmon modes. The right 
panels (5), (d) and (/) show the numerical solutions for the Landau damped (dashed blue lines) and undamped (red curves) 
plasmon branches. The plots show results for various distances between the surface and the layers: aikp = 2.0, 0.2 and 0.5, 
respectively. The second layer is placed at a distance a 2 , equal to twice as large as ai. All the plots are provided for extrinsic 
graphene (doped) with zero energy bandgap. 


and (q, uj) denotes the three dimensional (3D) bulk dielectric function of the thick-slab material. Furthermore, 
Eq. (IT^ introduces the definitions 


K: 


3D 


= i / 


COS qzZ 


0 


dqz 


(15) 


and 


cos qzZ 


2 

= - / dq^ - 2 , 2 ^ 3Dr -^ 

Making use of these results in Eq. ([II|), we then obtain 


(16) 
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FIG. 3: (Color online) Density plots showing the bulk modes in the miniband continuum and the line traces of plasmon dispersion 
for N = 4 graphene layers on a conducting substrate. The line plots show damped (dashed blue lines) and undamped (red 
curves) plasmon excitations. In {a)-{d), the layers are gapless, and in (e)-(/i) each layer has a gap A = 0.6 /i. The layers 
are equally spaced with inter-layer spacing ai 2 = . The separation between the first layer and the surface was chosen as 

di kp = 2.0 ai kp = ii/hvp. If the gap or number of layers is increased, the lowest branch does not re-appear for large q\\. It is 
more complicated than we thought. 


Gt/|| 


t(0) 


I + e 


-9ii —ihiil. 


1 + e(g||) 




1 + e{q\\) 


poo 

Jo 


(17) 


If the bulk plasma within the semi-infinite slab is fully local in the sense that e^(q, cj) eB{^) = 1 — in 

terms of the bulk plasma frequency cjp, then we use 




^oo = 


i(yi 

e-9iiN'l 

9||£b(w) 


(18) 


from which the corresponding local inverse dielectric function (qy, 2 ;, 0 '; w) may be obtained using Eq. (13) 


While Eq. (17) is useful because of its apparent simplicity, it is necessary to understand that its validity is restricted 
because the g^-integrations in Eqs. (p^ through extend over an infinite integration range. This blends the effects 

of the boundary/image length scale with that of the g;^-nonlocality dependence, eliminating the possibility of an 
unrestricted limit g^ ^ 0 and modifying the plasmon dependence on gy. Additionally, the imaginary part of (q, uj) 
is accounted for in these g;^-integrations, which consequently contributes to damping of the surface plasmon modes 
even in the low-gy limit. The ’’nonlocal” gy-correction to the surface plasmon, and its imaginary part involving 
damping, depend on the properties of the bounding surface. But, there is a range of applicability. This can be 
seen by examining the parameter measuring the importance of nonlocality in the bulk dielectric function 5^(q, cj), 

namely ^ {Tn"^/Ethermai) , where the characteristic thermal energy Ethermai = /i is the Eermi energy in the 
degenerate substrate with electron effective mass m*. Eor gy << it is reasonable to neglect nonlocality, at least 
in the surface plasmon frequency ujg = uJp/V2 as well as in comparison with other more pertinent sources of nonlocal 
behavior (but it cannot be neglected in the damping of the surface plasmon). 
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FIG. 4: (Color online) Plasmon dispersion for = 7 graphene layers on a conducting substrate. As in Fig. the layers in 
(a)-(d) are gapless, whereas in (e)-(/i) each layer has a gap A = 0.6 fi. The inter-layer separation is ai 2 = kp^. The hrst layer is 
at distances ai kp = 0.2 from the surface of the semi-infinite conductor. If the gap or number of layers is increased, the lowest 
branch does not re-appear for large q \\. 


III. NUMERICAL RESULTS 


In this section, we present numerical results for the plasmon dispersion for a system consisting of a semi-infinite 
conducting medium which is Coulomb coupled to N = 1,2,4 layers of graphene as shown in Figs, [^through We 
note that both the plasmon solutions and damping by bulk modes in the miniband continuum crucially depend on the 
separation between the constituents as well as the energy gap between the valence and conduction bands. For a single 
layer, our results shown in Fig. [^demonstrate that if the plasmon mode enters a region with Im n 2 ^(g||,co’) 7 ^ 0, the 
mode is Landau damped. Our calculations also show that when the distance a is less than a critical value dc — 0.4/c^^, 
in terms of the Fermi wave vector ii/hvp^ the lower acoustic plasmon mode is over-damped and this behavior seems 
analogous to data reported experimentally IMS]. This is obviously the case if the plasmon branch goes below the 
main diagonal uj = vf(1\\‘ The damping, as well as the critical distance changes in the presence of an energy bandgap 
for graphene. 

Similar conclusions for a pair of graphene layers electrostatically coupled to a semi-infinite conducting material are 
presented in Figj^ The principal difference between the case when there are two Coulomb-coupled layers is that if the 
distance of the layer nearest the conductor is less than the critical separation dc, both symmetric and antisymmetric 
modes become damped, for different ranges of wave vector. We emphasize that the upper plasmon branch (symmetric 
mode) remains almost unchanged for all cases, either with one or two graphene layers. 

The role played by the energy band gap is an important part of our investigation. For monolayer graphene, an 
energy gap leads to an extended region of undamped plasmons m- As we mentioned before, we pay particular 
attention to the regions outside of the single-particle excitation continuum with = 0 , since the plasmons 
























in these regions are not Landau-damped. In Figs. and we have plotted the plasmon dispersion relation for = 4 
and N = 7 graphene layers without and with an energy bandgap as well as for various distances between the nearest 
layer to the conducting surface. These results show that for a conducting substrate surface plasmon frequency denoted 
by Co’s = the surface plasmon frequency ujc of the hybrid superlattice always lies below Co’s. Furthermore, the 

intensity of this mode depends on the distance of the graphene layers from the surface of the conductor as well as the 
energy band gap between the valence and conduction bands of MLG. In the absence of a gap, our calculation shows 
that when the conductor surface-MLG separation exceeds a critical distance dc, the intensity of the surface plasmon in 
the long wavelength regime is high and may be detected up to some cut-off wave vector Qc. For the intensity 

of uJc is very weak until the plasmon wave vector exceeds some value q'^. However, when the surface-MLG separation 
is less than dc, the surface plasmon intensity in the long wavelength regime is weak and the mode only appears at 
shorter wavelength when q\\ > q'^. For gapped graphene, the surface plasmon frequency uJc is completely suppressed 
when the surface-layer separation is less than dc. 


IV. CONCLUDING REMARKS 

The appearance of a surface plasmon polariton for a multi-layer structure consisting of 2DELs was predicted in the 
1980’s by Giuliani and Quinn [2] as well as by Jain and Allen [3]. This surface mode is free from Landau damping or 
damping by bulk modes in the miniband continuum and lies the continuum of bulk modes. Additionally, this 

surface mode only exists above a critical wave vector gy which is determined by the layer spacing and the difference 
in the background dielectric constants of the layers and the surrounding medium. When gy, there is damping by the 
bulk plasmon modes. As a matter of fact, the dispersion equation for the layered superlattice structures investigated 
in urn is a special case which may be obtained from our more general Eq. 0 where we included the effects arising 
from a substrate. 

Very recently, in a series of experiments to determine the nature and behavior of plasmon excitations in graphene 
interacting with metallic substrates, Politano, et al. leHU] showed how the dispersion and intensity may be affected 
in a substantial way. The experimental results show that self-doped graphene supported by a metal substrate has two 
plasmon branches. There is an acoustic plasmon, with a linear dispersion, and a nonlinear plasmon. Both plasmon 
branches are similar in nature to those we presented in Eig. originating from the presence of a substrate. The 
present paper investigating the effects of a substrate on the plasmon excitations in a superlattice of graphene was 
stimulated by the experimental results on supported graphene. We are aware of the theoretical work on free-standing 
superlattice structures of MLGs [7], but the results there do not address the coupling to a metallic substrate which 
drastically affect the plasmon dispersion relation as may be observed from Eigs. [2 through 

The important conclusions of our work are as follows. We formulated and exploited a newly derived expression for 
plasmon dispersion in a superlattice of 2D layers which are Goulomb-coupled to a metallic substrate by taking into 
account the full nonlocality of the layers as well as the underlying conductor. We predict the existence of low-frequency 
nonlocal plasmon excitations at the vacuum-surface interface for various conditions of surface-layer separation as well 
as the energy gap. When the separation between the conducting surface and the nearest layer is less than some 
critical distance dc, the surface plasmon may not exist in the long wavelength limit. We obtain a surface plasmon at 
both intermediate as well as long wavelengths as this layer-surface separation is increased. Eor this hybrid structure, 
the surface plasmon frequency lies below the surface plasmon frequency for the semi-infinite substrate. Experimental 
verification of these simulated results may be achieved using high-resolution electron-energy-loss spectroscopy (EELS) 
m, for example. This paper was inspired by recent experimental work investigating the effect due to a metal on the 
collective plasmon mode of a single layer of graphene mm- We presented a new approach for generating a tunable 
surface plasmon using hybrid semiconductors. Additionally, our proposed approach based on hybrid semiconductors 
can be generalized to include other novel two-dimensional materials, such as hexagonal boron nitride, molybdenum 
disulfide and tungsten diselenide. 
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